# Questioning More: RT, Outward Facing Propaganda, and the Post-West World Order
# Erin Baggott Carter and Brett Carter
# This file replicates all tables and figures in the main text and online appendix, in the order that they appear

library(stargazer)
library(likert)
library(arm) 
library(effects) 
library(MASS)
library(ltm)
library(lattice)
library(psych)
library(car)
library(pastecs)
library(scales)
library(ggplot2)
library(arules)
library(plyr)
library(Rmisc)
library(GPArotation)
library(gdata)
library(MASS)
library(qpcR)
library(dplyr)
library(gtools)
library(Hmisc)
library(glm.predict)
library(quanteda)
library(readtext)
library(extrafont)
#font_import() 
loadfonts()
font_import(paths = c("/Library/Application Support/Apple/Fonts/iWork/Garamond.ttc", prompt = F))
library(RColorBrewer)
library(scales)
library(wordcloud)
library(tm)
library(mediation)

dir = '~/Dropbox/Autocratic External Propaganda/Production Materials/Replication/'
figdir = '~/Dropbox/Autocratic External Propaganda/Figures/' 
setwd(figdir)

######################################################
#                    load data                       #
######################################################

df = read.csv(paste(dir,'News Consumption and Public Opinion Survey_April 29, 2018_17.57.csv',sep=''),header=T,stringsAsFactors=F)
df$age = 2018-as.integer(as.vector(df$yob))
df$age[df$age>=100] <- NA 
df$Duration..in.seconds.=as.integer(as.vector(df$Duration..in.seconds.))
df = df[df$nationality=="American",];nrow(df)

df$drop <- 0
for(i in 1:nrow(df)){
  if(df$Duration..in.seconds.[i]<(4*60) & df$FL_13[i]!="Control-Nonewsarticle"){df$drop[i]<-1}
  if(df$Duration..in.seconds.[i]<(3*60) & df$FL_13[i]=="Control-Nonewsarticle"){df$drop[i]<-1}
}
dat <- df
df <- df[df$drop==0,];nrow(df)


######################################################
#                     data prep                      #
######################################################

df$party = NA
df$party[df$partisanship %in% c("Strongly Republican","Republican","Lean Republican")] = "Republican"
df$party[df$partisanship %in% c("Strongly Democrat","Democrat","Lean Democrat")] = "Democrat"
df$party[df$partisanship %in% c("Independent")] = "Independent"
df$party = as.factor(df$party)

# make partisanship continuous (1 = rep, 7 = dem)
df$partisanship[df$partisanship==""] = NA
df$partisanship[df$partisanship=="Prefer not to answer"] = NA
df$partisanship[df$partisanship=="Strongly Republican"] = 1
df$partisanship[df$partisanship=="Republican"] = 2
df$partisanship[df$partisanship=="Lean Republican"] = 3
df$partisanship[df$partisanship=="Independent"] = 4
df$partisanship[df$partisanship=="Lean Democrat"] = 5
df$partisanship[df$partisanship=="Democrat"] = 6
df$partisanship[df$partisanship=="Strongly Democrat"] = 7
df$partisanship = as.numeric(df$partisanship)

# make education, income continuous
df$education[df$education=="Less than a high school diploma"] = 0
df$education[df$education=="High school degree or equivalent (e.g. GED)"] = 1
df$education[df$education=="Some college, no degree"] = 2
df$education[df$education=="Associate degree (e.g. AA, AS)"] = 3
df$education[df$education=="Bachelor’s degree (e.g. BA, BS)"] = 4
df$education[df$education=="Master’s degree (e.g. MA, MS, MEd)"] = 5
df$education[df$education=="Professional degree (e.g. MD, DDS, DVM)"] = 6
df$education[df$education=="Doctorate (e.g. PhD, EdD)"] = 7
df$education = as.numeric(df$education)

# make income continuous
df$income[df$income=="Prefer not to answer"] = NA
df$income[df$income=="Less than $20,000"] = 0
df$income[df$income=="$20,000 to $34,999"] = 1
df$income[df$income=="$35,000 to $49,999"] =  2
df$income[df$income=="$50,000 to $74,999"] = 3
df$income[df$income=="$75,000 to $99,999"] = 4
df$income[df$income=="$100,000 to $149,999"] =5 
df$income[df$income=="$150,000 to $199,999"] = 6
df$income[df$income=="$200,000 to $249,000"] = 7
df$income[df$income=="Over $250,000"] = 8
df$income = as.numeric(df$income)

# make gender dichotomous
df$gender[df$gender=="Other (please specify):"] = NA
df = df[is.na(df$gender)==F,]
df$gender = factor(df$gender,c("Female","Male"))

# collapse into 2 treatments 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy(ForeignAgent)"] = 'FA'
df$FL_13[df$FL_13=="Control-Nonewsarticle"] = "CONTROL" 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy"] = "NON-FA"
mylevels = c("CONTROL","NON-FA","FA")
df$FL_13 = factor(df$FL_13,mylevels)

# low quality news share
unique(sort(df$newsSources))
tokens = c("breitbart","drudgeReport","infowars","intercept","alternet","redState","worldTruth","patribotics")
df$garbagenews=0
df$garbagenewsshare = 0
for(i in 1:nrow(df)){
  for(t in tokens){
    if(grepl(t,tolower(df$newsSources[i]))==T){
      df$garbagenews[i] = 1
    }
  }
  totalSources = length(strsplit(df$newsSources[i],",")[[1]])
  df$garbagenewsshare[i] = df$garbagenews[i]/totalSources
}

holder = df

# median respondent
median(df$age,na.rm=T)
table(df$gender)
median(df$education,na.rm=T)
median(df$income,na.rm=T)
median(df$partisanship,na.rm=T)

###########################################################################
#     Table 1: Effect of RT on Domestic Policy Views and Russia Views     #
###########################################################################

###### presidential approval
df = holder
table(df$approval)
df = df[df$approval!="",]
mylevels = c("Strongly disapprove","Somewhat disapprove","Neither approve nor disapprove","Somewhat approve","Strongly approve")
df$approval = factor(df$approval,levels=mylevels)
approval = polr(approval ~ FL_13 + age + education + income + gender + partisanship  + garbagenewsshare,data=df, Hess=TRUE); summary(approval)

###### economic trajectory
df = holder
table(df$trajectory)
df = df[df$trajectory!="",]
mylevels = c("Getting worse","Getting better")
df$trajectory = factor(df$trajectory,levels=mylevels)
trajectory = glm(trajectory ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(trajectory)

###### newsDiscernment
df = holder
table(df$newsDiscernment)
df = df[df$newsDiscernment!="",]
mylevels = c("Extremely difficult", "Moderately difficult", "Slightly difficult", "Neither easy nor difficult", "Slightly easy", "Moderately easy", "Extremely easy")
df$newsDiscernment = factor(df$newsDiscernment,levels=mylevels)
discernment = polr(newsDiscernment ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(discernment)

###### trust
df = holder
table(df$trust)
df = df[df$trust!="",]
mylevels = c("Strongly disagree", "Disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Agree", "Strongly agree")
df$trust = factor(df$trust,levels=mylevels)
trust = polr(trust ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(trust)

###### russia
df = holder
table(df$russia)
df = df[df$russia!="Don't know",]
mylevels = c("Very unfavorably", "Somewhat unfavorably", "Neither favorably nor unfavorably", "Somewhat favorably", "Very favorably")
df$russia = factor(df$russia,levels=mylevels)
russia = polr(russia ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(russia)

###### adversary
df = holder
table(df$russiaAdversary)
df = df[df$russiaAdversary!="",]
mylevels = c("Not much of a problem", "A serious problem but not an adversary", "An adversary")
df$russiaAdversary = factor(df$russiaAdversary,levels=mylevels)
adversary = polr(russiaAdversary ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(adversary)

###### production table
stargazer(approval,trajectory,discernment,trust,russia,adversary,font.size="small",title="Effect of RT on Domestic Policy Views and Russia Views",label="Table: Effect")


###########################################################################
#              Table 2: Effect of RT on Foreign Policy Views              #
###########################################################################

###### active
df = holder
df = df[df$active!="",];nrow(df)
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
active = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(active)

###### isolation
df = holder
table(df$isolation)
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
worldproblems = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(worldproblems)

###### allies
df = holder
table(df$interests)
mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
interests = glm(interests ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(interests)

##### engagement index: active, isolation, interests

df = holder
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
df$active = as.integer(df$active)

mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
df$isolation = as.integer(df$isolation)

mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
df$interests = as.integer(df$interests)

dataset = df[c('active','isolation','interests')]

# Create matrix: some tests will require it
data_matrix = cor(dataset, use = 'complete.obs')
round(data_matrix, 2) # Appendix Table 5: Correlation Matrix
cortest.bartlett(dataset)
KMO(data_matrix)
psych::alpha(dataset, check.keys=T)$total$std.alpha

pc1 = psych::principal(dataset, nfactors=length(dataset), rotate="none")
pc1

pc2 = psych::principal(dataset, nfactors=2, rotate = "varimax", scores = TRUE)
pc2
pc2$loadings

# Healthcheck
pc2$residual
pc2$fit
pc2$communality

dataset = cbind(dataset, pc2$scores)
summary(dataset$RC1, dataset$RC2)

newdat <- cbind(df, pc2$scores[,1])
m1 = lm(pc2$scores[,1] ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)

new1 <- pc2$scores[,1] + -1*min(pc2$scores[,1],na.rm=T)
for( i in 1:length(new1)){
  new1[i] = new1[i]  / max(new1,na.rm=T)
}

newdat <- cbind(df, new1)
m1 = lm(new1 ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)
g=m1

holderpca = df

##### production table
stargazer(active, worldproblems, interests, g,font.size="small",title="Effect of RT on Foreign Policy Views",label="Table: Effect")


###########################################################################
#              Figure 4: The Effect of RT on Global Engagement            #
###########################################################################

###### active -- 3 level
df = holder
df$active[df$active=="Strongly disagree"] = "Disagree"
df$active[df$active=="Somewhat disagree"] = "Disagree"
df$active[df$active=="Strongly agree"] = "Agree"
df$active[df$active=="Somewhat agree"] = "Agree"
mylevels = c("Disagree","Neither agree nor disagree","Agree")
df$active = factor(df$active,levels=mylevels)

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Should US be More Active Abroad?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsActive3.pdf" , device=pdf)

exp(-0.38630)
exp(-0.86148)

##### solve world problems (isolation)

df = holder
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,levels=mylevels)

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="US Does Right Amount to Solve World Problems?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="Effectsisolation.pdf" , device=pdf)

exp(-0.477081)
exp(-0.570775)

###### interests

df = holder
table(df$interests)
df$interests=as.character(as.vector(df$interests))
df$interests[df$interests=="In world affairs, the US should take into account the interests of its allies"] = 0 ###
df$interests[df$interests=="In world affairs, the US should follow its own national interests"] = 1

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = glm(as.numeric(interests) ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Follow National Interests over Allies' Interests",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsInterests.pdf" , device=pdf)

exp(0.214827)
exp(0.561351)

###########################################################################
#               Figure 1: RT Content from April 2018                      #
###########################################################################

df = read.csv(paste(dir,'RT_com_tweets_labeled_v3.csv',sep=''),header=T)
df = df[df$label!="",]
df$label = factor(df$label)

x.labels <- c("Criticizing\nDemocracies", "US Foreign\nPolicy", "US Domestic\nPolitics", "Criticizing\nIsrael", "Viral\nContent", "Championing the\nGlobal South", "RT Credibility", "Russian\nPolitics", "Miscellaneous", "Syria")
df$label = as.character(as.vector(df$label))
df$label[df$label=="Criticizing Israel"] = "Criticizing democracies"
df$label = as.factor(df$label)

#pdf("rtContent4Garamond.pdf", family="EB Garamond", width=4, height=4.5)
par(mar=c(bottom=12, left=4, top=3, right=0))
barplot(sort(table(df$label), decreasing=T),col=c(rep("darkgrey",1),rep("grey89",1),rep("darkgrey",2), rep("grey89",5)),las=2, ylab="Counts",ylim=c(0,800))
#dev.off()

###########################################################################
#                Figure 2: RT Propaganda Narratives                       #
###########################################################################

###### scott's keyness   

dat <- readtext(paste(dir,'RT_com_tweets_labeled_v3.csv',sep=''), text_field = "text",encoding="UTF8",verbosity=3)
corpus <- corpus(dat) 
corpus$superclass = corpus$label
corpus$superclass[corpus$superclass %in% c("US domestic politics","US foreign policy","Criticizing democracies")==F] = "Other"
badterms = c("_t","th","#","オ","_s","s","op-ed","�","rtuknews","says")

##### USDP   

subcorpus <- corpus_subset(corpus, superclass %in% c("US domestic politics","Other"))
subdfm <- dfm(subcorpus, groups = "superclass", 
              remove = c(badterms,stopwords("english")),
              remove_punct=T,remove_numbers=T,tolower=T,stem=F,remove_url=T,
              what="word1", 
              ngrams=1:3
              )

word <- c('israeli', 'american','america', 'russian','syrian')
tran <- c('israel',  'us',      'us',      'russia', 'syria')
featnames(dfm_select(subdfm, pattern = word))
subdfm <- dfm_replace(subdfm, pattern = word, replacement = tran)
result_keyness <- textstat_keyness(subdfm, target = "US domestic politics")

#pdf("quanteda_usdp_1panel_garamond.pdf", family="EB Garamond", width=4, height=4.5)
textplot_keyness(result_keyness,show_reference = F,show_legend=F)
#dev.off()


##### USFP   

target = "US foreign policy"
subcorpus <- corpus_subset(corpus, superclass %in% c(target,"Other"))
subdfm <- dfm(subcorpus, groups = "superclass", 
              remove = c(badterms,stopwords("english")),
              remove_punct=T,remove_numbers=T,tolower=T,stem=F,remove_url=T,
              what="word1"
              )

word <- c('israeli', 'american','america', 'russian','syrian')
tran <- c('israel',  'us',      'us',      'russia', 'syria')
featnames(dfm_select(subdfm, pattern = word))
subdfm <- dfm_replace(subdfm, pattern = word, replacement = tran)
result_keyness <- textstat_keyness(subdfm, target = target)

#pdf("quanteda_usfp_1panel_garamond.pdf", family="EB Garamond", width=4, height=4.5)
textplot_keyness(result_keyness,show_reference=F,show_legend=F)
#dev.off()


##### Criticizing Democracies   

target = "Criticizing democracies"
subcorpus <- corpus_subset(corpus, superclass %in% c(target,"Other"))
subdfm <- dfm(subcorpus, groups = "superclass", 
              remove = c(badterms,stopwords("english")),
              remove_punct=T,remove_numbers=T,tolower=T,stem=F,remove_url=T,
              what="word1"
              )

word <- c('israeli', 'american','america', 'russian','syrian','german', "germany's", 'french','ukrainian','skripals','british','europe','swedish','migrant')
tran <- c('israel',  'us',      'us',      'russia', 'syria', 'germany','germany',   'france','ukraine',  'skripal', 'uk',     'eu',    'sweden', 'migrants')
featnames(dfm_select(subdfm, pattern = word))
subdfm <- dfm_replace(subdfm, pattern = word, replacement = tran)
result_keyness <- textstat_keyness(subdfm, target = target)

#pdf("quanteda_critdem_1panel_garamond.pdf", family="EB Garamond", width=4, height=4.5)
textplot_keyness(result_keyness,show_reference=F,show_legend=F)
#dev.off()


###########################################################################
#     Appendix Figure 2: Most Distinctive Terms in RT Coverage            #
###########################################################################

# these data were prepared in python with Kessler's scattertext algorithm

##### most distinctive terms: criticizing democracies

dat = read.csv(paste(dir,'critScattertextData.csv',sep=''),stringsAsFactors=F);head(dat)
dat = dat[order(dat$term),]
dat$term = tolower(dat$term)
dat = dat[dat$term %in% c("st","fm","pm","one","two","amp","go","doesnt","ron","th","il","bn",
                          "rt","rtuknews","is","did","get","dont","rts","say")==F,]

dat$preScale <- dat$precount / dat$predocs		
dat$postScale <- dat$postcount / dat$postdocs	

hi <- 0.70
lo <- 1 - hi
max <- 1.0
pre <- dat[dat$x<lo & dat$y>hi & dat$y<max,]; nrow(pre)
post <- dat[dat$x>hi & dat$x<max & dat$y<lo,]; nrow(post)
both <- dat[dat$x>hi & dat$y>hi,]
neither <- dat[dat$x<lo & dat$y<lo,]

sample <- seq(from=1, to=100, by=1)
input.pre <- na.omit(pre[order(-pre$preScale), c("preScale")][sample])
names <- na.omit(pre[order(-pre$preScale), c("term")][sample])

layout(matrix(c(1,2), nrow=2), heights=c(0.3,3))
par(mar=c(0,0,0,0))
plot.new()
text(x=0.5, y=0.6, "Most Distinctive Terms:", cex=1.5)
text(x=0.5, y=0.15, "Criticizing Democracies", cex=1.5)
set.seed(1234)
wordcloud(scale=c(4,1),words=names, freq=input.pre, min.freq=0, max.words=200, random.order=F, rot.per=0, colors=brewer.pal(8, "Dark2"),main=NA)
#dev.print(file='critScattertext.pdf', device=pdf)

##### most distinctive terms: us foreign policy

dat = read.csv(paste(dir,'usfpScattertextData.csv',sep=''),stringsAsFactors=F);head(dat)
dat = dat[order(dat$term),]
dat$term = tolower(dat$term)
dat = dat[dat$term %in% c("st","fm","pm","one","two","amp","go","doesnt","ron","th","il","bn",
                          "rt","rtuknews","is","did","get","dont","rts","say")==F,]

dat$preScale <- dat$precount / dat$predocs		
dat$postScale <- dat$postcount / dat$postdocs	

hi <- 0.70
lo <- 1 - hi
max <- 1.0
pre <- dat[dat$x<lo & dat$y>hi & dat$y<max,]; nrow(pre)
post <- dat[dat$x>hi & dat$x<max & dat$y<lo,]; nrow(post)
both <- dat[dat$x>hi & dat$y>hi,]
neither <- dat[dat$x<lo & dat$y<lo,]

sample <- seq(from=1, to=100, by=1)
input.pre <- na.omit(pre[order(-pre$preScale), c("preScale")][sample])
names <- na.omit(pre[order(-pre$preScale), c("term")][sample])

layout(matrix(c(1,2), nrow=2), heights=c(0.3,3))
par(mar=c(0,0,0,0))
plot.new()
text(x=0.5, y=0.6, "Most Distinctive Terms:", cex=1.5)
text(x=0.5, y=0.15, "US Foreign Policy", cex=1.5)
set.seed(1234)
wordcloud(scale=c(4,1),words=names, freq=input.pre, min.freq=0, max.words=200, random.order=F, rot.per=0, colors=brewer.pal(8, "Dark2"),main=NA)
#dev.print(file='usfpScattertext.pdf', device=pdf)

##### most distinctive terms: us domestic politics 

dat = read.csv(paste(dir,'usdpScattertextData.csv',sep=''),stringsAsFactors=F);head(dat)
dat = dat[order(dat$term),]
dat$term = tolower(dat$term)
dat = dat[dat$term %in% c("st","fm","pm","one","two","amp","go","doesnt","ron","th","il","bn",
                          "rt","rtuknews","is","did","get","dont","rts","say","ad")==F,]

dat$preScale <- dat$precount / dat$predocs		
dat$postScale <- dat$postcount / dat$postdocs	

hi <- 0.70
lo <- 1 - hi
max <- 1.0
pre <- dat[dat$x<lo & dat$y>hi & dat$y<max,]; nrow(pre)
post <- dat[dat$x>hi & dat$x<max & dat$y<lo,]; nrow(post)
both <- dat[dat$x>hi & dat$y>hi,]
neither <- dat[dat$x<lo & dat$y<lo,]

sample <- seq(from=1, to=100, by=1)
input.pre <- na.omit(pre[order(-pre$preScale), c("preScale")][sample])
names <- na.omit(pre[order(-pre$preScale), c("term")][sample])

layout(matrix(c(1,2), nrow=2), heights=c(0.3,3))
par(mar=c(0,0,0,0))
plot.new()
text(x=0.5, y=0.6, "Most Distinctive Terms:", cex=1.5)
text(x=0.5, y=0.15, "US Domestic Politics", cex=1.5)
set.seed(1234)
wordcloud(scale=c(4,1),words=names, freq=input.pre, min.freq=0, max.words=200, random.order=F, rot.per=0, colors=brewer.pal(8, "Dark2"),main=NA)
#dev.print(file='usdpScattertext.pdf', device=pdf)

###########################################################################
# Appendix Figure 3: Demographic Characteristics of the Survey Population #
###########################################################################

df = read.csv(paste(dir,'News Consumption and Public Opinion Survey_April 29, 2018_17.57.csv',sep=''),header=T,stringsAsFactors=F)
df$age = 2018-as.integer(as.vector(df$yob))
df$age[df$age>=100] <- NA 
unique(df$age)
df$Duration..in.seconds.=as.integer(as.vector(df$Duration..in.seconds.))
df = df[df$nationality=="American",];nrow(df)

df$drop <- 0
for(i in 1:nrow(df)){
  if(df$Duration..in.seconds.[i]<(4*60) & df$FL_13[i]!="Control-Nonewsarticle"){df$drop[i]<-1}
  if(df$Duration..in.seconds.[i]<(3*60) & df$FL_13[i]=="Control-Nonewsarticle"){df$drop[i]<-1}
}
dat <- df
df <- df[df$drop==0,];nrow(df)


##### partisanship
mylevels = c("","Prefer not to answer","Strongly Republican","Republican","Lean Republican","Independent","Lean Democrat","Democrat","Strongly Democrat")
df$partisanship = factor(df$partisanship,levels=mylevels)
levels(df$partisanship) = mylevels = c("Don't know","Private","Strong R","Republican","Lean R","Independent","Lean D","Democrat","Strong D")
par(mar=c(14,5,4,2)+0.1)
barplot(table(df$partisanship),beside=T, space=0,axes=F,border=T,col.axis="White",xlab=NA,ylab=NA,main="Partisanship",col="grey",cex.main=2.5)
axis(1,at=c(1:9)-0.5,levels(df$partisanship),par(las=2),cex=2.5,cex.axis=2.5)
axis(2,cex=2.5,cex.axis=2.5)
#dev.print(file="BarplotPartisanship.pdf" , device=pdf)

##### age
df$age.rev[df$age<20] <- 1
df$age.rev[df$age>=20 & df$age<30] <- 2
df$age.rev[df$age>=30 & df$age<40] <- 3
df$age.rev[df$age>=40 & df$age<50] <- 4
df$age.rev[df$age>=50 & df$age<60] <- 5
df$age.rev[df$age>=60 & df$age<70] <- 6
df$age.rev[df$age>=70 & df$age<80] <- 7
df$age.rev[df$age>=80 & df$age<90] <- 8

mylevels = c("< 20",
             "20s",
             "30s",
             "40s",
             "50s", 
             "60s",
             "70s",
             "80s")

par(mar=c(14,5,4,2)+0.1)
barplot(table(df$age.rev),beside=T,space=0,main="Age",axes=F,col.axis="White",xlab=NA,ylab=NA,cex=2.5,cex.axis=2.5,cex.main=2.5)
axis(1,at=c(1:8)-0.5,mylevels,par(las=2),cex=2.5,cex.axis=2.5)
axis(2,cex=2.5,cex.axis=2.5)
#dev.print(file="BarplotAge.pdf" , device=pdf)

##### education
mylevels = c("Less than a high school diploma",
             "High school degree or equivalent (e.g. GED)",
             "Some college, no degree",
             "Associate degree (e.g. AA, AS)", 
             "Bachelor’s degree (e.g. BA, BS)",
             "Professional degree (e.g. MD, DDS, DVM)",
             "Master’s degree (e.g. MA, MS, MEd)",
             "Doctorate (e.g. PhD, EdD)")

df$education = factor(df$education,levels=mylevels)
levels(df$education) = c("< High school",
                         "High school",
                         "Some college",
                         "Associate", 
                         "Bachelor",
                         "Professional",
                         "Master",
                         "Doctorate")

par(mar=c(14,5,4,2)+0.1)
barplot(table(df$education),beside=T,space=0,main="Education",col.axis="White",axes=F,xlab=NA,ylab=NA,cex=2.5,cex.axis=2.5,cex.main=2.5)
axis(1,at=c(1:8)-0.5,levels(df$education),par(las=2),cex=2.5,cex.axis=2.5)
axis(2,cex=2.5,cex.axis=2.5)
#dev.print(file="BarplotEducation.pdf" , device=pdf)

##### income

mylevels = c("Prefer not to answer",
             "Less than $20,000",
             "$20,000 to $34,999",
             "$35,000 to $49,999", 
             "$50,000 to $74,999",
             "$75,000 to $99,999",
             "$100,000 to $149,999",
             "$150,000 to $199,999",
             "$200,000 to $249,000",
             "Over $250,000"
)

df$income = factor(df$income,levels=mylevels)
levels(df$income) = c("Private",
                      "<20",
                      "20-35",
                      "35-50", 
                      "50-75",
                      "75-100",
                      "100-150",
                      "150-200",
                      "200-250",
                      ">250"
)

par(mar=c(14,5,4,2)+0.1)
barplot(table(df$income),beside=T,space=0,main="Income ($1000)",cex=2.5,cex.axis=2.5,axes=F,xlab=NA,ylab=NA,cex.main=2.5)
axis(2,cex=2.5,cex.axis=2.5)
#dev.print(file="BarplotIncome.pdf" , device=pdf)


##### gender 

mylevels = c("Female","Male","Other")
df$gender = factor(df$gender,levels=mylevels)
levels(df$gender) = mylevels 

par(mar=c(14,5,4,2)+0.1)
barplot(table(df$gender),beside=T,space=0,main="Gender",axes=F,xlab=NA,ylab=NA,cex=2.5,cex.axis=2.5,cex.main=2.5)
axis(2,cex=2.5,cex.axis=2.5)
#dev.print(file="BarplotGender.pdf" , device=pdf)

#################################################################################################
#   Appendix Table 2: Effect of RT on Foreign Policy Views (Disaggregated Treatment Groups)     #
#################################################################################################

df = read.csv(paste(dir,'News Consumption and Public Opinion Survey_April 29, 2018_17.57.csv',sep=''),header=T,stringsAsFactors=F)
df$age = 2018-as.integer(as.vector(df$yob))
df$age[df$age>=100] <- NA
df$Duration..in.seconds.=as.integer(as.vector(df$Duration..in.seconds.))
df = df[df$nationality=="American",];nrow(df)

df$drop <- 0
for(i in 1:nrow(df)){
  if(df$Duration..in.seconds.[i]<(4*60) & df$FL_13[i]!="Control-Nonewsarticle"){df$drop[i]<-1}
  if(df$Duration..in.seconds.[i]<(3*60) & df$FL_13[i]=="Control-Nonewsarticle"){df$drop[i]<-1}
}
dat <- df
df <- df[df$drop==0,];nrow(df)


df$party = NA
df$party[df$partisanship %in% c("Strongly Republican","Republican","Lean Republican")] = "Republican"
df$party[df$partisanship %in% c("Strongly Democrat","Democrat","Lean Democrat")] = "Democrat"
df$party[df$partisanship %in% c("Independent")] = "Independent"
df$party = as.factor(df$party)

# for partisanship, make continuous (1 = rep, 7 = dem)
df$partisanship[df$partisanship==""] = NA
df$partisanship[df$partisanship=="Prefer not to answer"] = NA
df$partisanship[df$partisanship=="Strongly Republican"] = 1
df$partisanship[df$partisanship=="Republican"] = 2
df$partisanship[df$partisanship=="Lean Republican"] = 3
df$partisanship[df$partisanship=="Independent"] = 4
df$partisanship[df$partisanship=="Lean Democrat"] = 5
df$partisanship[df$partisanship=="Democrat"] = 6
df$partisanship[df$partisanship=="Strongly Democrat"] = 7
df$partisanship = as.numeric(df$partisanship)

# make education, income continuous
df$education[df$education=="Less than a high school diploma"] = 0
df$education[df$education=="High school degree or equivalent (e.g. GED)"] = 1
df$education[df$education=="Some college, no degree"] = 2
df$education[df$education=="Associate degree (e.g. AA, AS)"] = 3
df$education[df$education=="Bachelor’s degree (e.g. BA, BS)"] = 4
df$education[df$education=="Master’s degree (e.g. MA, MS, MEd)"] = 5
df$education[df$education=="Professional degree (e.g. MD, DDS, DVM)"] = 6
df$education[df$education=="Doctorate (e.g. PhD, EdD)"] = 7
df$education = as.numeric(df$education)

df$income[df$income=="Prefer not to answer"] = NA
df$income[df$income=="Less than $20,000"] = 0
df$income[df$income=="$20,000 to $34,999"] = 1
df$income[df$income=="$35,000 to $49,999"] =  2
df$income[df$income=="$50,000 to $74,999"] = 3
df$income[df$income=="$75,000 to $99,999"] = 4
df$income[df$income=="$100,000 to $149,999"] =5 
df$income[df$income=="$150,000 to $199,999"] = 6
df$income[df$income=="$200,000 to $249,000"] = 7
df$income[df$income=="Over $250,000"] = 8
df$income = as.numeric(df$income)

# make gender dichotomous
df$gender[df$gender=="Other (please specify):"] = NA
df = df[is.na(df$gender)==F,]
df$gender = factor(df$gender,c("Female","Male"))

# all treatments
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy(ForeignAgent)"] = "DEM-FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics(ForeignAgent)"] = "USDP-FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy(ForeignAgent)"] = 'USFP-FA'
df$FL_13[df$FL_13=="Control-Nonewsarticle"] = "CONTROL" 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy"] = "DEM"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics"] = "USDP"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy"] = "USFP"
mylevels = c("CONTROL","DEM","DEM-FA","USDP","USDP-FA","USFP","USFP-FA")
df$FL_13 = factor(df$FL_13,mylevels)

# check low quality news
unique(sort(df$newsSources))
tokens = c("breitbart","drudgeReport","infowars","intercept","alternet","redState","worldTruth","patribotics")
df$garbagenews=0
df$garbagenewsshare = 0
for(i in 1:nrow(df)){
  for(t in tokens){
    if(grepl(t,tolower(df$newsSources[i]))==T){
      df$garbagenews[i] = 1
    }
  }
  totalSources = length(strsplit(df$newsSources[i],",")[[1]])
  df$garbagenewsshare[i] = df$garbagenews[i]/totalSources
}


holder = df

###### active
df = holder
df = df[df$active!="",];nrow(df)
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
active = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(active)

###### isolation
df = holder
table(df$isolation)
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
worldproblems = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(worldproblems)

###### allies
df = holder
table(df$interests)
mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
interests = glm(interests ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(interests)

##### engagement index: active, isolation, interests

df = holder
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
df$active = as.integer(df$active)

mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
df$isolation = as.integer(df$isolation)

mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
df$interests = as.integer(df$interests)

dataset = df[c('active','isolation','interests')]

# Create matrix: some tests will require it
data_matrix = cor(dataset, use = 'complete.obs')
round(data_matrix, 2)
cortest.bartlett(dataset)
KMO(data_matrix)
psych::alpha(dataset, check.keys=T)$total$std.alpha

pc1 = psych::principal(dataset, nfactors=length(dataset), rotate="none")
pc1

pc2 = psych::principal(dataset, nfactors=2, rotate = "varimax", scores = TRUE)
pc2
pc2$loadings

# Healthcheck
pc2$residual
pc2$fit
pc2$communality

dataset = cbind(dataset, pc2$scores)
summary(dataset$RC1, dataset$RC2)

newdat <- cbind(df, pc2$scores[,1])
m1 = lm(pc2$scores[,1] ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)

new1 <- pc2$scores[,1] + -1*min(pc2$scores[,1],na.rm=T)
for( i in 1:length(new1)){
  new1[i] = new1[i]  / max(new1,na.rm=T)
}

newdat <- cbind(df, new1)
m1 = lm(new1 ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)
g=m1

# back up pca
holderpca = origholderpca = df

##### production table
stargazer(active, worldproblems, interests, g,font.size="small",title="Effect of RT on Foreign Policy Views (Disaggregated Treatment Groups)",label="Table: Effect")


############################################################################################
# Appendix Figure 4: The Effect of RT on Global Engagement, Disaggregated Treatment Groups #
############################################################################################

##### active -- 3 level
df = holder
df$active[df$active=="Strongly disagree"] = "Disagree"
df$active[df$active=="Somewhat disagree"] = "Disagree"
df$active[df$active=="Strongly agree"] = "Agree"
df$active[df$active=="Somewhat agree"] = "Agree"
mylevels = c("Disagree","Neither agree nor disagree","Agree")
df$active = factor(df$active,levels=mylevels)

m1 = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Should the US be More Active Abroad?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsActive7.pdf" , device=pdf)


##### solve world problems (isolation)

df = holder
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,levels=mylevels)

m1 = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Does US do Right Amount to Solve World Problems?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="Effectsisolation_7.pdf" , device=pdf)

###### interests

df = holder
table(df$interests)
df$interests=as.character(as.vector(df$interests))
df$interests[df$interests=="In world affairs, the US should take into account the interests of its allies"] = 0 ###
df$interests[df$interests=="In world affairs, the US should follow its own national interests"] = 1

m1 = glm(as.numeric(interests) ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Follow National Interests over Allies' Interests",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsInterests7.pdf" , device=pdf)

#################################################################################################
#   Appendix Table 3: Effect of RT on Foreign Policy Views (Disaggregated by Party Affiliation) #
#################################################################################################

df = read.csv(paste(dir,'News Consumption and Public Opinion Survey_April 29, 2018_17.57.csv',sep=''),header=T,stringsAsFactors=F)
df$age = 2018-as.integer(as.vector(df$yob))
df$age[df$age>=100] <- NA  
df$Duration..in.seconds.=as.integer(as.vector(df$Duration..in.seconds.))
df = df[df$nationality=="American",];nrow(df)

df$drop <- 0
for(i in 1:nrow(df)){
  if(df$Duration..in.seconds.[i]<(4*60) & df$FL_13[i]!="Control-Nonewsarticle"){df$drop[i]<-1}
  if(df$Duration..in.seconds.[i]<(3*60) & df$FL_13[i]=="Control-Nonewsarticle"){df$drop[i]<-1}
}
dat <- df
df <- df[df$drop==0,];nrow(df)


df$party = NA
df$party[df$partisanship %in% c("Strongly Republican","Republican","Lean Republican")] = "Republican"
df$party[df$partisanship %in% c("Strongly Democrat","Democrat","Lean Democrat")] = "Democrat"
df$party[df$partisanship %in% c("Independent")] = "Independent"
df$party = as.factor(df$party)

# for partisanship, make continuous (1 = rep, 7 = dem)
df$partisanship[df$partisanship==""] = NA
df$partisanship[df$partisanship=="Prefer not to answer"] = NA
df$partisanship[df$partisanship=="Strongly Republican"] = 1
df$partisanship[df$partisanship=="Republican"] = 2
df$partisanship[df$partisanship=="Lean Republican"] = 3
df$partisanship[df$partisanship=="Independent"] = 4
df$partisanship[df$partisanship=="Lean Democrat"] = 5
df$partisanship[df$partisanship=="Democrat"] = 6
df$partisanship[df$partisanship=="Strongly Democrat"] = 7
df$partisanship = as.numeric(df$partisanship)

# make education, income continuous
df$education[df$education=="Less than a high school diploma"] = 0
df$education[df$education=="High school degree or equivalent (e.g. GED)"] = 1
df$education[df$education=="Some college, no degree"] = 2
df$education[df$education=="Associate degree (e.g. AA, AS)"] = 3
df$education[df$education=="Bachelor’s degree (e.g. BA, BS)"] = 4
df$education[df$education=="Master’s degree (e.g. MA, MS, MEd)"] = 5
df$education[df$education=="Professional degree (e.g. MD, DDS, DVM)"] = 6
df$education[df$education=="Doctorate (e.g. PhD, EdD)"] = 7
df$education = as.numeric(df$education)

df$income[df$income=="Prefer not to answer"] = NA
df$income[df$income=="Less than $20,000"] = 0
df$income[df$income=="$20,000 to $34,999"] = 1
df$income[df$income=="$35,000 to $49,999"] =  2
df$income[df$income=="$50,000 to $74,999"] = 3
df$income[df$income=="$75,000 to $99,999"] = 4
df$income[df$income=="$100,000 to $149,999"] =5 
df$income[df$income=="$150,000 to $199,999"] = 6
df$income[df$income=="$200,000 to $249,000"] = 7
df$income[df$income=="Over $250,000"] = 8
df$income = as.numeric(df$income)

# make gender dichotomous
df$gender[df$gender=="Other (please specify):"] = NA
df = df[is.na(df$gender)==F,]
df$gender = factor(df$gender,c("Female","Male"))

# collapse into 2 treatments 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy(ForeignAgent)"] = 'FA'
df$FL_13[df$FL_13=="Control-Nonewsarticle"] = "CONTROL" 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy"] = "NON-FA"
mylevels = c("CONTROL","NON-FA","FA")
df$FL_13 = factor(df$FL_13,mylevels)

# check low quality news
unique(sort(df$newsSources))
tokens = c("breitbart","drudgeReport","infowars","intercept","alternet","redState","worldTruth","patribotics")
df$garbagenews=0
df$garbagenewsshare = 0
for(i in 1:nrow(df)){
  for(t in tokens){
    if(grepl(t,tolower(df$newsSources[i]))==T){
      df$garbagenews[i] = 1
    }
  }
  totalSources = length(strsplit(df$newsSources[i],",")[[1]])
  df$garbagenewsshare[i] = df$garbagenews[i]/totalSources
}

holder = df

###### active
df = holder
df = df[df$active!="",];nrow(df)
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
active = polr(active ~ FL_13 + age + education + income + gender + party + garbagenewsshare,data=df, Hess=TRUE); summary(active)

###### isolation
df = holder
table(df$isolation)
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
worldproblems = polr(isolation ~ FL_13 + age + education + income + gender + party + garbagenewsshare,data=df, Hess=TRUE); summary(worldproblems)

###### allies
df = holder
table(df$interests)
mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
interests = glm(interests ~ FL_13  + age + education + income + gender + party + garbagenewsshare,data=df, family="binomial"); summary(interests)

##### engagement index: active, isolation, interests

df = holder
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
df$active = as.integer(df$active)

mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
df$isolation = as.integer(df$isolation)

mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
df$interests = as.integer(df$interests)

dataset = df[c('active','isolation','interests')]

data_matrix = cor(dataset, use = 'complete.obs')
round(data_matrix, 2)
cortest.bartlett(dataset)
KMO(data_matrix)
psych::alpha(dataset, check.keys=T)$total$std.alpha

pc1 = psych::principal(dataset, nfactors=length(dataset), rotate="none")
pc1

pc2 = psych::principal(dataset, nfactors=2, rotate = "varimax", scores = TRUE)
pc2
pc2$loadings

# Healthcheck
pc2$residual
pc2$fit
pc2$communality

dataset = cbind(dataset, pc2$scores)
summary(dataset$RC1, dataset$RC2)

newdat <- cbind(df, pc2$scores[,1])
m1 = lm(pc2$scores[,1] ~ FL_13  + age + education + income + gender + party + garbagenewsshare, data=newdat); summary(m1)

new1 <- pc2$scores[,1] + -1*min(pc2$scores[,1],na.rm=T) 
for( i in 1:length(new1)){
  new1[i] = new1[i]  / max(new1,na.rm=T)
}

newdat <- cbind(df, new1)
m1 = lm(new1 ~ FL_13  + age + education + income + gender + party + garbagenewsshare, data=newdat); summary(m1)
g=m1

# back up pca
holderpca = df

##### production table
stargazer(active, worldproblems, interests, g,font.size="small",title="Effect of RT on Foreign Policy Views (Disaggregated Treatment Groups)",label="Table: Effect")

############################################################################################
#       Appendix Figure 5: The Effect of RT on Global Engagement, by Party Membership      #
############################################################################################

# active -- 3 level
df = holder
df$active[df$active=="Strongly disagree"] = "Disagree"
df$active[df$active=="Somewhat disagree"] = "Disagree"
df$active[df$active=="Strongly agree"] = "Agree"
df$active[df$active=="Somewhat agree"] = "Agree"
mylevels = c("Disagree","Neither agree nor disagree","Agree")
df$active = factor(df$active,levels=mylevels)

m1 = polr(active ~ FL_13 + age + education + income + gender + party + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13","party"), mod = m1), rug = FALSE,style="stacked",main="Should the US be More Active Abroad?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsActiveParty.pdf" , device=pdf)

#### solve world problems (isolation)

df = holder
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,levels=mylevels)

m1 = polr(isolation ~ FL_13 + age + education + income + gender + party + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13","party"), mod = m1), rug = FALSE,style="stacked",main="Does US do Right Amount to Solve World Problems?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsIsolationParty.pdf" , device=pdf)

###### interests

df = holder
table(df$interests)
df$interests=as.character(as.vector(df$interests))
df$interests[df$interests=="In world affairs, the US should take into account the interests of its allies"] = 0 ###
df$interests[df$interests=="In world affairs, the US should follow its own national interests"] = 1

m1 = glm(as.numeric(interests) ~ FL_13 + age + education + income + gender + party + garbagenewsshare,data=df, family="binomial"); summary(m1)
plot(Effect(focal.predictors = c("FL_13","party"), mod = m1), rug = FALSE,style="stacked",main="Follow National Interests over Allies' Interests",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsInterestsParty.pdf" , device=pdf)

######################################################################################################
# Appendix Table 4: Effect of RT on Foreign Policy Views (5 Minutes Treatment, 4 Minutes Control)    #
######################################################################################################

df = read.csv(paste(dir,'News Consumption and Public Opinion Survey_April 29, 2018_17.57.csv',sep=''),header=T,stringsAsFactors=F)
df$age = 2018-as.integer(as.vector(df$yob))
df$age[df$age>=100] <- NA 
df$Duration..in.seconds.=as.integer(as.vector(df$Duration..in.seconds.))
df = df[df$nationality=="American",];nrow(df)

df$drop <- 0
for(i in 1:nrow(df)){
  if(df$Duration..in.seconds.[i]<(5*60) & df$FL_13[i]!="Control-Nonewsarticle"){df$drop[i]<-1}
  if(df$Duration..in.seconds.[i]<(4*60) & df$FL_13[i]=="Control-Nonewsarticle"){df$drop[i]<-1}
}
dat <- df
df <- df[df$drop==0,];nrow(df)


df$party = NA
df$party[df$partisanship %in% c("Strongly Republican","Republican","Lean Republican")] = "Republican"
df$party[df$partisanship %in% c("Strongly Democrat","Democrat","Lean Democrat")] = "Democrat"
df$party[df$partisanship %in% c("Independent")] = "Independent"
df$party = as.factor(df$party)

# for partisanship, make continuous (1 = rep, 7 = dem)
df$partisanship[df$partisanship==""] = NA
df$partisanship[df$partisanship=="Prefer not to answer"] = NA
df$partisanship[df$partisanship=="Strongly Republican"] = 1
df$partisanship[df$partisanship=="Republican"] = 2
df$partisanship[df$partisanship=="Lean Republican"] = 3
df$partisanship[df$partisanship=="Independent"] = 4
df$partisanship[df$partisanship=="Lean Democrat"] = 5
df$partisanship[df$partisanship=="Democrat"] = 6
df$partisanship[df$partisanship=="Strongly Democrat"] = 7
df$partisanship = as.numeric(df$partisanship)

# make education, income continuous
df$education[df$education=="Less than a high school diploma"] = 0
df$education[df$education=="High school degree or equivalent (e.g. GED)"] = 1
df$education[df$education=="Some college, no degree"] = 2
df$education[df$education=="Associate degree (e.g. AA, AS)"] = 3
df$education[df$education=="Bachelor’s degree (e.g. BA, BS)"] = 4
df$education[df$education=="Master’s degree (e.g. MA, MS, MEd)"] = 5
df$education[df$education=="Professional degree (e.g. MD, DDS, DVM)"] = 6
df$education[df$education=="Doctorate (e.g. PhD, EdD)"] = 7
df$education = as.numeric(df$education)

df$income[df$income=="Prefer not to answer"] = NA
df$income[df$income=="Less than $20,000"] = 0
df$income[df$income=="$20,000 to $34,999"] = 1
df$income[df$income=="$35,000 to $49,999"] =  2
df$income[df$income=="$50,000 to $74,999"] = 3
df$income[df$income=="$75,000 to $99,999"] = 4
df$income[df$income=="$100,000 to $149,999"] =5 
df$income[df$income=="$150,000 to $199,999"] = 6
df$income[df$income=="$200,000 to $249,000"] = 7
df$income[df$income=="Over $250,000"] = 8
df$income = as.numeric(df$income)

# make gender dichotomous
df$gender[df$gender=="Other (please specify):"] = NA
df = df[is.na(df$gender)==F,]
df$gender = factor(df$gender,c("Female","Male"))

# collapse into 2 treatments 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics(ForeignAgent)"] = "FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy(ForeignAgent)"] = 'FA'
df$FL_13[df$FL_13=="Control-Nonewsarticle"] = "CONTROL" 
df$FL_13[df$FL_13=="Treatment-CriticizingDemocracy"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USDomesticPolitics"] = "NON-FA"
df$FL_13[df$FL_13=="Treatment-USForeignPolicy"] = "NON-FA"
mylevels = c("CONTROL","NON-FA","FA")
df$FL_13 = factor(df$FL_13,mylevels)

# check low quality news
unique(sort(df$newsSources))
tokens = c("breitbart","drudgeReport","infowars","intercept","alternet","redState","worldTruth","patribotics")
df$garbagenews=0
df$garbagenewsshare = 0
for(i in 1:nrow(df)){
  for(t in tokens){
    if(grepl(t,tolower(df$newsSources[i]))==T){
      df$garbagenews[i] = 1
    }
  }
  totalSources = length(strsplit(df$newsSources[i],",")[[1]])
  df$garbagenewsshare[i] = df$garbagenews[i]/totalSources
}

holder = df

###### active
df = holder
df = df[df$active!="",];nrow(df)
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
active = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(active)

###### isolation
df = holder
table(df$isolation)
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
worldproblems = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(worldproblems)

###### allies
df = holder
table(df$interests)
mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
interests = glm(interests ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(interests)

##### engagement index: active, isolation, interests

df = holder
mylevels = c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Agree","Strongly agree")
df$active = factor(df$active,mylevels)
df$active = as.integer(df$active)

mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,mylevels)
df$isolation = as.integer(df$isolation)

mylevels = c("In world affairs, the US should take into account the interests of its allies","In world affairs, the US should follow its own national interests")
df$interests = factor(df$interests,mylevels)
df$interests = as.integer(df$interests)

dataset = df[c('active','isolation','interests')]

data_matrix = cor(dataset, use = 'complete.obs')
round(data_matrix, 2)
cortest.bartlett(dataset)
KMO(data_matrix)
psych::alpha(dataset, check.keys=T)$total$std.alpha

pc1 = psych::principal(dataset, nfactors=length(dataset), rotate="none")
pc1

pc2 = psych::principal(dataset, nfactors=2, rotate = "varimax", scores = TRUE)
pc2
pc2$loadings

# Healthcheck
pc2$residual
pc2$fit
pc2$communality

dataset = cbind(dataset, pc2$scores)
summary(dataset$RC1, dataset$RC2)

newdat <- cbind(df, pc2$scores[,1])
m1 = lm(pc2$scores[,1] ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)

new1 <- pc2$scores[,1] + -1*min(pc2$scores[,1],na.rm=T) 
for( i in 1:length(new1)){
  new1[i] = new1[i]  / max(new1,na.rm=T)
}

newdat <- cbind(df, new1)
m1 = lm(new1 ~ FL_13  + age + education + income + gender + partisanship + garbagenewsshare, data=newdat); summary(m1)
g=m1

# back up pca
holderpca = df

##### production table
stargazer(active, worldproblems, interests, g,font.size="small",title="Effect of RT on Foreign Policy Views (Disaggregated Treatment Groups)",label="Table: Effect")


##############################################################################################################
#   Appendix Figure 6: The Effect of RT on Global Engagement (More Restrictive Completion Time Threshold)    #
##############################################################################################################

# active -- 3 level
df = holder
df$active[df$active=="Strongly disagree"] = "Disagree"
df$active[df$active=="Somewhat disagree"] = "Disagree"
df$active[df$active=="Strongly agree"] = "Agree"
df$active[df$active=="Somewhat agree"] = "Agree"
mylevels = c("Disagree","Neither agree nor disagree","Agree")
df$active = factor(df$active,levels=mylevels)

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = polr(active ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Should US be More Active Abroad?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsActive3Restricted.pdf" , device=pdf)

#### solve world problems (isolation 4)

df = holder
df = df[df$isolation!="",];nrow(df)
mylevels = c("Does too much","Does right amount","Does too little")
df$isolation = factor(df$isolation,levels=mylevels)

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = polr(isolation ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, Hess=TRUE); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="US Does Right Amount to Solve World Problems?",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsisolationRestricted.pdf" , device=pdf)

###### interests

df = holder
table(df$interests)
df$interests=as.character(as.vector(df$interests))
df$interests[df$interests=="In world affairs, the US should take into account the interests of its allies"] = 0 ###
df$interests[df$interests=="In world affairs, the US should follow its own national interests"] = 1

df$FL_13 = as.character(as.vector(df$FL_13))
df$FL_13[df$FL_13=="CONTROL"] = "Control"
df$FL_13[df$FL_13=="NON-FA"] = "Nondisclosure"
df$FL_13[df$FL_13=="FA"] = "Disclosure"
mylevels = c("Control","Nondisclosure","Disclosure")
df$FL_13 = factor(df$FL_13,levels=mylevels)

m1 = glm(as.numeric(interests) ~ FL_13 + age + education + income + gender + partisanship + garbagenewsshare,data=df, family="binomial"); summary(m1)
plot(Effect(focal.predictors = c("FL_13"), mod = m1), rug = FALSE,style="stacked",main="Follow National Interests over Allies' Interests",ylab="Response Probability",xlab="Treatment Condition")
#dev.print(file="EffectsInterestsRestricted.pdf" , device=pdf)

###########################################################################
#              Appendix Table 6: Causal Mediation Analysis                #
###########################################################################

df = origholderpca

##### trajectory
df.med <- na.omit(newdat[newdat$FL_13=="CONTROL" | newdat$FL_13=="NON-FA" | newdat$FL_13=="FA",c("new1", "FL_13", "trajectory","partisanship", "age", "education", "income", "gender","garbagenewsshare")])
df.med$treatment <- 0
df.med$treatment[df.med$FL_13!="CONTROL"] <- 1
df.med$trajectoryRev <- 0
df.med = df.med[df.med$trajectory!="",]
df.med$trajectoryRev = 0
df.med$trajectoryRev[df.med$trajectory=="Getting worse"] = 1
med.fit = lm(trajectoryRev ~ treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(med.fit)
out.fit = lm(new1 ~ trajectoryRev + treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(out.fit)
med.out <- mediate(med.fit, out.fit, treat = "treatment", mediator = "trajectoryRev", sims = 100, boot=TRUE)
summary(med.out)
a=med.fit
b=out.fit

##### approval
df.med <- na.omit(newdat[newdat$FL_13=="CONTROL" | newdat$FL_13=="NON-FA" | newdat$FL_13=="FA",c("new1", "FL_13", "approval","partisanship", "age", "education", "income", "gender","garbagenewsshare")])
df.med$treatment <- 0
df.med$treatment[df.med$FL_13!="CONTROL"] <- 1
df.med$approvalRev <- 0
df.med = df.med[df.med$approval!="",]
mylevels = c("Strongly disapprove","Somewhat disapprove","Neither approve nor disapprove","Somewhat approve","Strongly approve")
df.med$approval = factor(df.med$approval,levels=mylevels)
df.med$approvalRev = as.integer(df.med$approval)
med.fit = lm(approvalRev ~ treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(med.fit)
out.fit = lm(new1 ~ approvalRev + treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(out.fit)
med.out <- mediate(med.fit, out.fit, treat = "treatment", mediator = "approvalRev", sims = 100, boot=TRUE)
summary(med.out)
c=med.fit
d=out.fit

##### trust
df.med <- na.omit(newdat[newdat$FL_13=="CONTROL" | newdat$FL_13=="NON-FA" | newdat$FL_13=="FA",c("new1", "FL_13", "trust","partisanship", "age", "education", "income", "gender","garbagenewsshare")])
df.med$treatment <- 0
df.med$treatment[df.med$FL_13!="CONTROL"] <- 1
df.med$trustRev <- 0
df.med = df.med[df.med$trust!="",]
mylevels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree")
df.med$trust = factor(df.med$trust,levels=mylevels)
df.med$trustRev = as.integer(df.med$trust)
med.fit = lm(trustRev ~ treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(med.fit)
out.fit = lm(new1 ~ trustRev + treatment + age + education + income + gender + partisanship + garbagenewsshare, data=df.med); summary(out.fit)
med.out <- mediate(med.fit, out.fit, treat = "treatment", mediator = "trustRev", sims = 10000, boot=TRUE)
summary(med.out)
e=med.fit
f=out.fit

stargazer(a,b,c,d,e,f,font.size="small",title="Causal Mediation Analysis",label="Table: CMA")







